#!/usr/bin/perl
# 
# Exercise 11.5
# 
# Write a program that extracts the secondary structure information contained in the 
# HELIX, SHEET, and TURN record types of PDB files. Print out the secondary 
# structure and the primary sequence together, so that it's easy to see by what sec-
# ondary structure a given residue is included. (Consider using a special alphabet for 
# secondary structure, so that every residue in a helix is represented by H, for exam-
# ple.)
#
# Answer to Exercise 11.5
# This exercise involves writing subroutines to parse the HELIX, SHEET, and TURN lines.
# The main complication involves keeping track of chain names for these lines and for
# the SEQRES lines: we also have to modify the parseSEQRES subroutine for this.
# Note that if there is only one chain, the chain name is often left blank.

use strict;
use warnings;
use BeginPerlBioinfo;

# Give the name of your PDB file
my $pdbfile = 'pdb/44/pdb244l.ent';

# Get the file data, parse the record types
my @pdbfiledata = get_file_data($pdbfile);

my %recordtypes = parsePDBrecordtypes(@pdbfiledata);

# Extract the name and the sequence of each chain
my %chains = extractSEQRES2($recordtypes{'SEQRES'});
foreach my $key (keys %chains) {
	# Change the sequence from 3-letter to 1-letter amino acid codes
	$chains{$key} = iub3to1($chains{$key});
	# For debugging purposes
	#print "$key:$chains{$key}\n";
}

# Parse the HELIX, SHEET, and TURN record types
# The following three statement use the ternary operator CONDITION : IFTRUE : IFFALSE
# They test if the hash is defined for the key, and if so, parse it; if not, they set
# it to the empty list
$recordtypes{'HELIX'} ? my %helix = extractHELIX($recordtypes{'HELIX'}) : ();
$recordtypes{'SHEET'} ? my %sheet = extractSHEET($recordtypes{'SHEET'}) : ();
$recordtypes{'TURN'}  ? my %turn  = extractTURN ($recordtypes{'TURN'})  : ();

# For debugging purposes:
#foreach my $key (sort keys %helix) {
#	print "$key:$helix{$key}\n";
#}
#foreach my $key (sort keys %sheet) {
#	print "$key:$sheet{$key}\n";
#}
#foreach my $key (sort keys %turn) {
#	print "$key:$turn{$key}\n";
#}

# 
# Now make a annotation strings that contain the helix, sheet, and turn information
#

my %annotation = ();

foreach my $chain_name (sort keys %chains) {

	# Make a blank annotation string the same length as the sequence of the chain
	$annotation{$chain_name} = ' ' x length($chains{$chain_name});

	if( defined $helix{$chain_name} ) {
		foreach my $structure ( split /:/, $helix{$chain_name} ) {
			my($structure, $position) = split /;/, $structure;
			substr( $annotation{$chain_name}, $position - 1, length($structure)) = $structure;
		}
	}
	if( defined $sheet{$chain_name} ) {
		foreach my $structure ( split /:/, $sheet{$chain_name} ) {
			my($structure, $position) = split /;/, $structure;
			substr( $annotation{$chain_name}, $position - 1, length($structure)) = $structure;
		}
	}
	if( defined $turn{$chain_name} ) {
		foreach my $structure ( split /:/, $turn{$chain_name} ) {
			my($structure, $position) = split /;/, $structure;
			substr( $annotation{$chain_name}, $position - 1, length($structure)) = $structure;
		}
	}

	print "Chain name = \"$chain_name\"\n";
	print $chains{$chain_name},"\n";
	print "$annotation{$chain_name}\n";
}

exit;

################################################################################
# Subroutines
################################################################################

# extractHELIX
#
#-given an scalar containing HELIX lines,
#    return a scalar containing H

sub extractHELIX {
    my($helix) = @_;

    # make array of lines
    my @record = split ( /\n/, $helix);

    # hash to store chains
    my %chainhash = ();
    
    foreach my $line (@record) {
        # Chain is in column 20, starting position in columns 22-25
	#  length is in columns 72-76
	# "stripspaces" removes leading and trailing spaces
	# If there's only one chain in the PDB entry, the chain name may be blank
        my($thischain)	= stripspaces(substr($line, 19, 1));
        my($start)	= stripspaces(substr($line, 21, 4));
        my($length)	= stripspaces(substr($line, 71, 5));
    
	if( defined $chainhash{$thischain} ) {
	    $chainhash{$thischain} .= ':' . 'H' x $length . ";$start";
	}else{
	    $chainhash{$thischain} = 'H' x $length . ";$start";
	}
    }

    return %chainhash;
}

# extractSHEET
#
#-given an scalar containing SHEET lines,
#    return a scalar containing H

sub extractSHEET {
    my($sheet) = @_;

    # make array of lines
    my @record = split ( /\n/, $sheet);

    # hash to store chains
    my %chainhash = ();
    
    foreach my $line (@record) {
        # Chain is in column 22, start position is in columns 23-26
	#  end position is in columns 34-37
	# "stripspaces" removes leading and trailing spaces
        my($thischain)	= stripspaces(substr($line, 21, 1));
        my($start)	= stripspaces(substr($line, 22, 4));
        my($end)	= stripspaces(substr($line, 33, 4));
	my($length)	= $end - $start + 1;
    
	if( defined $chainhash{$thischain} ) {
	    $chainhash{$thischain} .= ':' . 'S' x $length . ";$start";
	}else{
	    $chainhash{$thischain} = 'S' x $length . ";$start";
	}
    }

    return %chainhash;
}

# extractTURN
#
#-given an scalar containing TURN lines,
#    return a scalar containing H

sub extractTURN {
    my($turn) = @_;

    # make array of lines
    my @record = split ( /\n/, $turn);

    # hash to store chains
    my %chainhash = ();
    
    foreach my $line (@record) {
        # Chain is in column 20, starting position in columns 21-24
	#  length is in columns 32-35
	# "stripspaces" removes leading and trailing spaces
        my($thischain)	= stripspaces(substr($line, 19, 1));
        my($start)	= stripspaces(substr($line, 20, 4));
        my($end)	= stripspaces(substr($line, 31, 4));
	my($length)	= $end - $start + 1;
    
	if( defined $chainhash{$thischain} ) {
	    $chainhash{$thischain} .= ':' . 'T' x $length . ";$start";
	}else{
	    $chainhash{$thischain} = 'T' x $length . ";$start";
	}
    }

    return %chainhash;
}


sub stripspaces {
	my($string) = @_;

	$string =~ s/^\s*//;
	$string =~ s/\s*$//;

	return $string;
}


# extractSEQRES2
#
#-given an scalar containing SEQRES lines,
#    return an array containing the chains of the sequence
#
# Modified version that reports the chain name plus the sequence, 
# returning the values as a hash

sub extractSEQRES2 {

    use strict;
    use warnings;

    my($seqres) = @_;

    my $lastchain;
    my $sequence = '';
    my %results = (  );
    # make array of lines

    my @record = split ( /\n/, $seqres);
    
    foreach my $line (@record) {
        # Chain is in column 12, residues start in column 20
        my ($thischain) = stripspaces(substr($line, 11, 1));
        my($residues)  = substr($line, 19, 52);
    
        # Check if a new chain, or continuation of previous chain
        if(not defined $lastchain ) {
            $sequence = $residues;
        }elsif("$thischain" eq "$lastchain") {
            $sequence .= $residues;
        # Finish gathering previous chain (unless first record)
        }elsif ( $sequence ) {
            $results{$lastchain} = $sequence;
            $sequence = $residues;
        }
        $lastchain = $thischain;
    }

    # save last chain
    $results{$lastchain} = $sequence;
    
    return %results;
}
